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q ! ABSTRACT 

' The dynamical evolution of binaries of intermediate-massive black holes (IMBHs, 

massive black holes with a mass ranging between 10 2 and 10 4 M Q ) in stellar clusters 
has recently received an increasing amount of attention. This is at least partially due 
to the fact that if the binary is hard enough to evolve to the phase at which it will 
start emitting gravitational waves (GWs) efficiently, there is a good probability that 
I/-) ' it will be detectable by future space-borne detectors like LISA. We study this evolu- 

, tion in the presence of rotation in the cluster by carrying out a series of simulations 

f*"*. ■ of an equal-mass binary of IMBHs embedded in a stellar distribution with different 

' rotational parameters. The survey indicates that eccentricities and inclinations are 

qq , primarily determined by the initial conditions of the IMBHs and the influence of dy- 

■ namical friction, even though they are finally perturbed by the scattering of field stars. 
Q\ | In particular, the eccentricity is strongly connected to the initial IMBHs velocities, and 
c"^) . values of ~ 0.7 up to 0.9 are reached for low initial velocities, while almost circular 

orbits result if the initial velocities are increased. Evidence suggests a dependency of 
. ^ | the eccentricity on the rotation parameter. We found only weak changes in the in- 

i clination, with slight variations of the orientation of the angular momentum vector 

' of the binary. Counter-rotation simulations yield remarkably different results in ec- 

centricity. A Monte Carlo study indicates that these sources will be detectable by a 
detector such as LISA with median signal to noise ratios of between 10 and 20 over a 
three year period, although some events had signal to noise ratios of 300 or greater. 
Furthermore, one should also be able to estimate the chirp-mass with median frac- 
tional errors of 10~ 4 , reduced mass on the order of 10~ 3 and luminosity distance on 
the order of 10 . Finally, these sources will have a median angular resolution in the 
LISA detector of about 3 square degrees, putting events firmly in the field of view of 
future electromagnetic detectors such as LSST. 

Key words: stellar dynamics - black hole physics - gravitational waves - detection 
- globular clusters: general 



1 INTRODUCTION 



* E-mail: Pau.Amaro-Seoane@aei.mpg.de 

f E-mail: Eichliorn@irs.uni-stuttgart.dc 

f E-mail: porter@apc.univ-paris7.fr 

§ E-mail: Spurzem@ari.uni-heidelberg.de 



Even though their existence are not as well-established as 
the existence of stellar-mass or supermassive black holes 
(SMBHs), it is plausible that intermediate-mass black holes 
(IMBHs; masses M ~ 10 2 ~ 4 M B ) exist i n the centre of 
stellar clusters fsee iMiller fc Colbertl [2004. and references 
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therein). The coalescence of a binary of two massive black 
holes in that mass range is a powerful source of gravity 
waves which will be detectable by space-born observatories 
such as the Laser Interferometer Space Antenna (LISAn) 
l|Amaro-Seoane fc Freitadl2006l ; lAmaro-Seoane et aLll2009n . 

There are two possible ways to theoretically ex- 
plain the presence of a binary of IMBHs in a globu- 
lar clu ster. (1) The "single cluster channel" : iGiirkan et al.l 
(2006) address this possibility in the scenario of a run- 
away growth of two stars in a young cluster via phys- 
ical collisions in the innermost part of the stellar sys- 
tem, whe re the heaviest stars have sunk through mass seg- 
regat i on dPorteeies Zwart fc McMillan! l2000l: IGiirkan et al.l 
120041 ; IPortegies Zwart et al.l 120041 : iFreitag et al.ll2006l ). By 
adding a fra c tion o f primordial binaries to the cluster, 
IGiirkan et al.l (|2006l ) found not one but two independent 
very massive stars growing in the centre. Eventually, and 
due to post-Newtonian inestabilities, they could collapse 
and form a massive black hole in the relevant mass range. 
lAmaro-Seoane et all (|2009l ) followed the evolution of such 
a binary with direct-summation N— body simulations and 
estimated that, in some cases, one can expect a resid- 
ual eccentricity as high as 0.3 for the binary when it en- 
ters the LISA band. (2 ) The "double cluster channel" : 
lAmaro-Seoane fc Freitad HH) discussed an alternative 
way to form binaries of IMBHs in stellar clusters. It has been 
observed that in star forming regions such as the Antennae or 
Stephan's Quintet, hundreds of young massive star clusters 
are clustered into larger complexes of a few 100 pc across 
dWhitmore et al.lll999l : IZhang fc Fallll999l : iGallagher et~ai] 
2001). These clusters contain typically some ~ 10 s stars 
within a fe w parsecs and it is m ost likely that some of them 
are bound ( Di eball et al.l [2002). Also, most of the clusters in 
binaries are coeval and younger than 300 Myr, which m eans 
that they merge early. lAmaro- Seoane fc Freitad (|2006h fol- 
lowed the evolution of two IMBHs initially located at the 
centre of two such clusters which collide, and studied the 
orbital decay to the centre. 

The evolution of the orbital parameters of the binary of 
massive BHs is determined by the stellar dynamics and the 
emission of GWs, which is negligible at long distances but 
becomes more and more dominant as the semi-major axis of 
the binary shrinks. In order to understand the distribution 
of parameters we can expect for these systems when they 
enter the window of detection of LISA, it is important to 
analyse the dynamical story of the binary in detail. The 
evolution starting at distances at which dynamical friction 
is important down to the phase of strong emission of GWs 
-and their detection and characterisation- is a complex and 
long process which requires different techniques to address 
it. 

One can distinguish roughly three different regimes in 
the process: (i) at the beginning, the massive BHs are at 
distances in which GWs, though always present, are totally 
negligible and the evolution is dominated purely by the dy- 
namics. At this stage, dynamical friction will sink the two 
massive BHs down to the centre. The perturbers, the sin- 
gle massive BHs in our case, are moving through a sea of 
small stars (as compared to their masses). The velocity vec- 



tor of the stars are rotated after deflecting with the IMBHs. 
The projected component in the direction of the deflection 
is shorter. Hence, the massive object, the IMBH, is cumulat- 
ing just after it a high-density stellar region. The perturber 
will feel a drag from that region from the conservation of 
J in the direction of its velocity vector. The direction does 
not change to first-order, but the amplitude decreases A drag 
force starts to act on to the perturber, so that it slows downs 
as it sinks down to the centre of the stellar system. This force 
happens to be proportional the square of the mass of the per- 
turber so, the bigger the mass of the perturber, the bigger 
the dynamical effects, in spite of the bigger inertia, (ii) As 
they approach closer and closer, the two massive BHs form 
a bound state, a binary system. At this stage, the binary 
interacts strongly with stars coming from the surrounding 
stellar system. Since the stars have a much smaller mass, 
the outcome of the interaction is that a star is slingshot into 
the stellar system with a higher kinetic energy, gained from 
the removal of energy and angular momentum of the binary; 
therefore, the semi-major axis of the binary, a, shrinks a bit 
more. These interactions mostly tend to increase the eccen- 
tricity of the binary (iii) The process is repeated again and 
again, provided the reservoir of stars is not empty - in which 
case the loss-cone would be depleted -, until the separation 
between the members of the binary is small enough that 
the emission of GWs is strong enough as to take over the 
dynamics as the main factor of evolution of the orbital pa- 
rameters. The binary is practically isolated from the stellar 
system. Thereafter, the binary begins to circularise. Obvi- 
ously, the transition between these phases in the evolution is 
not well-separated and th e whole evolution requires numer - 
ical tools to investigate it. Amaro-Seoane fc Freitad (2006); 
lAmaro-Seoane et al.l (|200E ) addressed some of these ques- 
tions in the context of globular clusters and IMBHs. They 
proved, with A^— body models, that slingshot ejections of 
stars increase the eccentricity of the IMBH binary to ~ 0.8 
and beyond and that later the emission of GWs then cir- 
cularises the orbit to rather low, yet detectable values of 
eccentricity. 

But Nature is more complex than that. A key effect 
that will determine the global evolution of both the cluster 
and of the massive binary of IMBHs is the rotation of the 
stellar system. Rotation has been identified in clusters for 
a long time. Deviations from sph erical symmetr y were dis- 
covered in some globular clusters (|Shaplevll 19301 ) as early as 
the beginning of the last century. This flattening is a finger- 
print for rotation and the measurements of ellipticity were 
later ex tended to galactic and extragalactic globul ar clusters 
(see e.g. lWhite fc Shawllll987l ; IStaneva et alj|l996h . One can 
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also detect rotation by measuring r adial velocities of indi- 
vidual stars in the globular clusters ijMevlan fc Mavor|[l986l ; 
iGebhardt et al.ll20 0ol; iReiins et al.ll2006l ). ~ 

lAmaro-Seoane fc Freitagl 1 20061 ) first addressed this 
problem focusing on the detection of GWs in the context 
of a binary of IMB H formed as the res ult of a collision of 
two clusters; lAmaro-Seoane et al.l (|2009l ) looked at the same 
problem from the perspective of a born-in binary, extended 
the study to multi-mass clusters and non-equal mass bina- 
ries and described the global dynamical evolution with the 
implication of the larger eccentricities they found in a gen- 
eral context of detection. 
area/indexThfml-fiieealdfeteff eccentricity of the binary is decisive 
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in the study. Direct summation JV-body models, including 
post-Newtonian corrections, show that the stellar dynam- 
ical history before the relativistic regime can significantly 
affect the final evolution an d leads to different merger 
times (|Berentzen et al.l [2009), In particular, it turns out 
that massive binaries may enter the relativistic phase with 
high eccentricities, and signatures of the eccentricities are 
kept in the harmonics of th e gravitational waveforms until 
the moment of coale s cence llAmaro-Seoane fc Freitaell2006l ; 



lAmaro-Seoane et alj 120091 ; iBerentzen et al.l 120091 ). The 



evolution of the eccentrici ty has been previously discussed 
in a n umber of articles dMakino et a] 1 19931 ; iHemsendorl 
iMnosavlievic fc MerrittJ l200ll; I? 



2000 



2002 



Berczik et al.l 



Berentzen et al 



2005; 



2009; [Berczik et all 120061 ; 



iHemsendorf et a] 



Amaro-Scoane & Frcitag 
lAmaro-Seoane et al.l [2009 ; 

(|l964T ). the orbit-averaged 
timescale of coalescence due to the emission of gravitational 
radiation is given by 



20091 ). 

As derived by iPetersI 



tew = 



c 5 4w 



64 G 3 M 1 M 2 (M 1 + M 2 )F{e) 



(1) 



where Mi, M2 denote the black hole masses, acw the char- 
acteristic separation for gravitational wave emission, G the 
gravitational constant, c the speed of light and 



F(e) = (l-e 2 ) 



-7/2 



1 73 2 
1+ 24 6 



37 4 
96 6 



(2) 



is a function with dependence on the eccentricity e. Thus the 
coalescence time can shrink by several orders of magnitude 
if the eccentricity is high enough, resulting in a stronger 
burst of gravitational radiation and different characteristic 
waveforms. Furthermore, the behavior of the the inclination 
of the orbital plane is potentially interesting in predicting 
processes related to angular momentum exchange between 
IMBHs and field stars. 

A natural continuation of the analysis carried out until 
now is to add another physical factor to the problem, the 
rotation of the system, since it can have a very important 
impact in the global dynamics of the cluster, ft can also 
particularly effect the evolution of the eccentricity of the 
binary and, thus, the detection and characterisation of the 
GW observation. 

From a standpoint of the data analysis of such GWs 
for LISA, because of their low mass -as compared to SMBH 
binaries-, IMBH binaries should be visible at moderate to 
high frequencies in the LISA detector. As most of these bina- 
ries coalesce outside the LISA band, they will be observable 
in the detector throughout the lifetime of the mission. This 
should allow us to confidently detect and estimate the pa- 
rameters for such systems. If such sources exist in the LISA 
data stream, it will also allow us, assuming a strong enough 
signal, to provide accurate distance measurements in the lo- 
cal universe. 

The stucture of this paper is as follows: We start by 
giving a description of the numerical method used for the 
simulations in section ([2} ; later, in section ([3]) we give a short 
overview of the initial conditions we use for the numerical 
study; in section @ we provide a detailed analysis of the dy- 
namics of the system, i.e. evolution of the binding energy, ec- 



centricity and inclination of the IMBH binary; in section (0) 
we study some cases in which the binary is initially set up 
on a counter-rotating orbit in relation to the stellar system 
in which it is embedded and study the associated Brownian 
motion; in section (0 we dicuss about the implications for 
lower-frequency Astrophysics and the detectability of such 
systems. In the last section ([8} we summarise the results and 
give the conclusions of the study. 



2 NUMERICAL METHOD 

The simulations in this work have been performed using 
NBODY6++, a . parallelised version of Aarseth's NBODY6 
jSpurzem|[l999T ). The code includes a Hermite integration 
scheme, KS-regularisation l|Kustaanheimp fc Stiefel 1965I) 
and t he Ahmad-Cohen neighbour scheme ( Ahmad fc Cohen] 
Il973h . No softening has been introduced; this circumstance 
allows an accurate treatment of the effects due to super- 
elastic scattering events, which play a crucial part in black 
hole binary evolution and require a precise calculation of the 
trajectories throughout the interaction. 

Additionally, in order to improve the exactness of com- 
putation of the motion of particles in the environment of 
the black holes, a modification in the determination of the 
neighbour radius in the Ahmad-Cohen scheme was imple- 
mented. In principle, the Ahmad-Cohen scheme divides the 
force on a particle into a regular and an irregular component, 
assigning both forces different time steps; the irregular com- 
ponent is computed over the nearest particles populating an 
area called the neighbour sphere. Normally, the neighbour 
sphere, which is in its extension defined by the neighbour ra- 
dius, is characterised as containing a defined number of stars 
rin, which is typically set to n n = 50 in simulations dealing 
with particle numbers of the order presented here. However, 
when considering a scenario consisting of two heavy parti- 
cles of the same mass, embedded in a stellar component of 
equal-mass particles, the neighbour radius is enlarged by a 
factor 



7 = 



- 1 



+ 1 



(3) 



if during the declaration of the neighbour particles of a par- 
ticle j and the ne ighbour candidat es i a mass difference 
mi/rrij 7^ 1 occurs (IHems endorl 200(J). The enlargement fac- 
tor 7 is symmetric in the masses rrn and m,j. As a result, 
massive particles (black holes) receive a larger neighbour ra- 
dius, and a massive particle is also more likely declared as 
neighbour of a stellar particle. This method accommodates 
the influence of a black hole on its surroundings and lessens 
the underestimation on that stellar component which pos- 
sessed a black hole nearby, but outside the neighbour sphere 
in the scheme without the enlargement factor. The param- 
eters j3 and A have been set /3 = 0.03 and A = 1, yield- 
ing 7 = 10.57, consisten t with the simulations presented by 
IHemsendorf et all l|2002l ). 
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0.3 


V2v c 
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0.3 


Vc 


0.256 


0.553 
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0.3 


0.136i> c 


0.256 


0.553 
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6 


0.6 


Vc 


0.320 


0.550 


19.81 




Q 


6 


0.6 


0.136t; c 


0.320 


0.550 


19.81 



Table 1. Overview over the initial conditions of all simulations. 
All values are expressed in N-body units. Wo is the King pa- 
rameter, uq the rotation parameter, vq the initial velocity of the 
IMBHs, r c the core radius, Vc the central velocity and Trot /^kin 
is the ratio of rotational energy and total kinetic energy at time 
T = 0. We use 7 = 73, (3 = 0.03, A = 1.0. In all simulations 
Vrcg = 0.02, whilst ?7in-cg = 0.02 for the cases A — T and 0.01 for 
the cases J — Q 



3 INITIAL CONDITIONS 

In all our simulations, the conversion factors are as follows: A 
unit of length is equivalent to U\^ B = 1.1 pc, a unit of mass 

= 13.233 
The set 



to U\l 



43921. 6Af©, a unit of velocity to U\{ 



km/s and a unit of time to W|m =8.14- 1CT 2 Myrs. 
of simulations was carried out for a total particle number 
N = 64000, including two massive black holes Mi = M2 = 
0.01 embedded in a dense stellar system of 63998 equal- 
mass particles m* = 1.5625 • 10 -5 WjM B - The initial stellar 
distribution satisfies a rotating King model, 



f(E,L z ) = 



const x e 



-/3(£-4>t) 





E < 



(4) 

where j3 = 1/cr 2 represents the inverse one-dimensional ve- 
locity dispersion, E and L z the energy and the z-component 
of the angular momentum of a star per unit mass, and 
$t the tidal potential of the model. The King parameters 
Wo and ojq are defined by Wo = — /3($o — $t) and ujq — 
^/9/A-KGpc ■ fi c with the parameters $0 representing the 
central potential, p c a mass densit y and Q c approximate! 
the angular velocity in the centre (lEinsel fc Spurzeml 199' 
iLaeoute fc Longarettilll996l : lLongaretti fc Lagoutdll996t l. 

A symmetric set-up underlies all simulations, with the 
rotation axis of the King model in the z-direction and the 
two black holes located in the xy-plane on the opposite side 
of the cluster on the core radius r c of the model. Thereby, 
the usual definition of the core radius in N-body simulations, 



N* 



2=1 



3=1 



1/2 



iV* ^ AT/5, (5) 



is used l|Aarsethl 120031 ') , where va = Y^j 
is the density centre of the system consisting of particles 
of the mass m with the coordinates rj, and pj = 3(fc — 
1)771/47 ^^^^ the local density in an area around each par- 
ticle j l|Casertano fc Hutlll985l ). The quantity r^ k ) is to be 
interpreted as the radius of the sphere over which the lo- 
cal density is evaluated, characterised in size by the number 
of stars k populating this volume. Applying this formalism, 
fc = 6 is the optimal choice. The summation in equation 
[5] is restricted to the innermost N/5 particles, which saves 
computation time while showing coevally adequate agree- 
ment with an exact calculations. The initial velocities are 
adjusted in the xy-plane tangential to a circle around the 
centre, adopting the values vo = 0.136i> c ,tj c and \^2v c in 
terms of the circular velocity v c . Scenarios with the black 
holes moving initially with, as well as contrary to, the rota- 
tion of the stellar system have been investigated. 

The evolution of the black hole binary was followed us- 
ing King parameters Wo = 3; 6 and wo = 0.0; 0.3; 0.6. We 
have chosen these parameters because they are representa- 
tive for the problem we want to address in this work, from 
the absence of rotation to a rather high value. Typical val- 
ues for wo in globular clusters range between 0.1 and 0.5. 
For instance, cuCen has a value of 0.5, N280 8 of 0.3, 47Tuc 
of 0.15, N5286 of 0.5 etc jFiestas et all2006h . 

The accuracy was tuned in such a way that relative en- 
ergy errors measured over one NBODY time were < 10~ 3 
concerning Wo = 3 and < 10 -4 concerning Wo = 6 simula- 
tions respectively. The complete survey of investigations is 
shown in Tab[T] 



4 DYNAMICS OF THE SYSTEM 

4.1 Evolution of the binding energy 

The total energy of the binary in a two-body approximation 
is given by 



GM1M2 



2pr 2 



(6) 



where r denotes their separation, Mi and M2 the black hole 
masses, [i — M\M-zf(M\ + M2) the reduced mass, I the 
angular momentum and G the gravitational constant. 

Fig. [T] shows the time evolution of the total energy. In 
both cases, the initial velocity is vo = 0.136u c on the core ra- 
dius, different colours represent different rotational param- 
eters. Naturally, as long as the gravitational force of the 
stellar system dominates the motion of the black holes, the 
two-body energy is not very meaningful. Initial oscillations 
are the result of this invalidity: Due to the symmetric set-up 
both black holes reach the apoapsis and the periapsis almost 
simultaneously, in the apoapsis (where their separation is at 
maximum, which consequently means a local minimum in 
the total energy) the black holes are formally bound to each 
other (E < 0). However, they feel the potential of the stellar 
system not included in eq[6]and are accelerated to the centre 
while gathering kinetic energy in such a way that the bound 
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Wo =3, Vo = 0.136 v c 



l(w = 0.6) 





Figure 1. Left panel: Evolution of the total energy of the IBMH binary for the models X, T and C (from the top to the bottom at later 
times; in the on-line version of the paper displayed in blue, green and red respectively) Right panel: As in the previous panel, for models 
Q, O and C 
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Figure 2. Left panel: Same as Fig. JTJ for models T>, £ and T. In this case all models have the same King parameter Wo = 3 and 
wo = 0.3, so that they differ only in the initial velocities (from the top to the bottom at later times vo = V2v c , vq = v c and vo = 0.136u c , 
depicted in the on-line version in red, green and blue, respectively) Right panel: As in the previous panel but with Wo = 6, for models 
M, Af and O 
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a i, 


H 
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0.526 


41.06 


210 


7.8±2.0 


T 


0.522 


43.53 


209 


7.2±1.9 


J 


0.478 


117.12 


191 


5.2±1.3 


K 


0.478 


117.12 


191 


5.6±1.5 


C 


0.478 


117.12 


191 


6.1±1.6 


N 


0.468 


108.04 


187 


5.8±1.5 


O 


0.468 


108.04 


187 


6.8±1.7 


Q 


0.425 


64.79 


170 


10.4±2.7 



Table 2. The hardening constants determined in simulations for 
a variety of models and initial conditions, uo represents the initial 
central velocity dispersion, no the initial central particle density, 
the characteristic separation for hardening via super-elastic 
scattering processes and H the hardening constant. 

state is resolved. In this first stage, each black hole individ- 
ually suffers dynamical friction, which is the main process 
of losing energy. 

The role of dynamical friction decreases when a per- 
manently bound state occurs (the energy remains negative), 
as the dynamical friction force acts primarily on the mo- 
tion of the now formed binary rather than on the individual 
black holes. Super-elastic scattering events of field stars at 
the binary become more and more important for the reduc- 
tion of its energy. These events become visible in the tiny- 



peak structure that appears in each curve in Fig[T]at times 
t > 17 for Wo = 3 and t > 8 for W = 6 respectively. In the 
stage where super-elastic scattering dominates the picture, 
the energy loss rate is commonly written in terms of the 
dimensionless hardening constant H 

±(1\=hg£ (7) 

where a is the separation of the black holes, p the mass 
density and a the velocity dispe rsion in the environment of 
the binary (see e.g. lMerrittil200ll ). The constant slope of the 
energy in Fig[T]is expected from eq[7]with p/a =const. 

In Tab(2] hardening constants have been determined for 
scenarios in which a constant energy loss rate had developed 
before the simulation ended. The time derivative of the in- 
verse separation was taken from the slope of the curves in 
Fig[T] which is connected to the energy by E — —Gp,/2a. 
Approximately, a/p — ao/po was assumed with po and uq 
as initial values within the 1% Lagrangian radius of the 
model. As this represents a rather vague approximation, a 
25% range of this ratio was combined with an regression 
estimate of the uncertainties in slope designation to obtain 
the error margins. 

Since the stage which is dominated by superelastic scat- 
tering is reached sooner if the central potential of the stellar 
distribution is deeper, the hardening constants have been 
calculated primarily for runs with Wo = 6 in Table 2. For 
these runs, the separation of the IMBHs has fallen below 
the characteristic separation au = GM/Aao at the end of 
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the simulation, which indicates that the system is in the 
hardening regime. 

The values of the hardening constants are in major- 
ity slightly below comp ared to the H = 8.4 published by 
iHemsendorf et all (|2002l ). where a Plummer model was used. 
The lower values can be possibly explained by the fact that 
dynamical friction might still have a noticeable influence. 
Regarding our calculated an as criteria for the domination 
of super-elastic scattering events, the hardening separation 
could be significantly smaller if a increases during the sim- 
ulation as ah oc <j~ 2 . An enhanced a can be expected for 
p/a =const. if it is assumed that the black hole would cap- 
ture stars during the simulation and raise the central density. 

We can see in Fig. 1 that King potentials of Wo = 3 as 
well as Wo — 6 the two-body energy reaches higher values 
(i.e. the system is less bound) at a given time for uo = 0.6, 
as compared to the simulations with u>o = 0.0 and luq = 0.3. 
Additionally, we find that in the presence of faster rotation 
loo = 0.6, the transfer of angular momentum to the field stars 
is inhibited during the first 2-4 time units. This can lead to 
the circularising of the IMBHs trajectories and, indeed, we 
find remarkable smaller eccentricities for this case (see Fig 
6 of next section). 

Variation of the initial velocities is presented in Fig[2] 
In simulations with wo = v c , no oscillations occur in the first 
time units, as the black holes spiral to the centre symmetri- 
cally on circular-like orbits. 



4.2 Eccentricity 

The eccentricity is determined by 



e= 4/1 



2EP 



FigED shows the complete survey over calculated eccen- 
tricity evolution. Each plot shows simulations of a fixed pair 
of King parameters under variation of the initial velocity. 
In the case of Wo — 6, luq = 0.6, no simulation could be 
performed with vq = %/2i> c without overstepping the error 
limits mentioned in section [3] holding the set of NBODY 
accuracy parameters. The runs were stopped when a fixed 
physical calculation time of the PC cluster was exceeded. 

Initial oscillations appear for the same reasons as in 
the plots of the total energy previously discussed. After the 
binary has been formed, and in principle represents a two 
body system perturbed by encountering field stars, the ec- 
centricity converges to a fixed value that underlies at most a 
weak drift. This behaviour is consistent wit h previous work 
jHemsendorfj |2000| ; IHemsendorf et"a3 12002| ). When the ec- 
centricity has swung into a certain level, again a tiny-peak 
structure develops as the result of super-elastic scattering 
processes of field stars at the hardening binary. Note that 
following the swing-in-procedure, the stochastic fluctuations 
are of the same order, due to the logarithmic scaling. 

All simulations displayed in Figure [3] with an initial ve- 
locity comparable to the circular velocity, tend to end up in 
low-eccentricity motions of the black hole components, while 
vo — 0.136v c runs reach generally higher fi nal eccentrici- 
ties. T his behaviour was already indicated bv lMakino et all 
{ 19931 ). who simulated two black holes of the masses M = 
0.01 in a Plummer sphere of 16348 particles. They found 



very high final eccentricities e ~ 0.99 applying very low ini- 
tial velocities, while their largest value, vo = 0.5v c , reached 
a noti ceably smaller fina l e ~ 0.665. Investigat i ons c arried 
out bv lHemsendorj (|2000l ) and IHemsendorf "e t al ( 20021) used 
initial velocities no = 0.136u c , their high eccentricities were 
verified in the simulations presented here for the rotating 
King model. 

The dependency of the final eccentricity on initial veloc- 
ities can be understood by considering the black hole trajec- 
tories. In Fig. 21 for vo = v c , the black holes spiral, at first 
independently of each other, to the centre. The influence 
of dynamical friction causes a steady loss of kinetic energy. 
Within the time interval t — [10.11; 20.14], the total energy 
becomes negative and the binary reaches a bound state; sub- 
sequently the binary hardens, the separation decreases due 
to super-elastic scattering events and the circular motion of 
the centre of mass of the binary itself becomes visible. At the 
time the attractive force between the black holes becomes 
comparable to the gravitational force of the stellar distribu- 
tion, the individual trajectories of the black holes are still 
circular around the systems centre of mass. This means that 
the circular orbits generated by the initial velocity is "con- 
served" until the binary reaches a bound state and beyond, 
since dynamical friction is not strong enough to change the 
trajectories dramatically. 

A different situation is obtained for vo = 0.136 v c . As a 
consequence of the low velocity, the black holes must plunge 
near to the centre, but dynamical friction is, at the time of 
the closest encounter (the periapsis of the relative motion), 
not sufficiently effective to prevent the re-swing to the outer 
regions and to circularise the orbits in this way. Therefore, 
the initial form of the orbits is kept until the end of the sim- 
ulation. The initial velocity Vo = V^v c is also non-circular. 
The deviations from the previous case are due to the fact 
that dynamical friction is stronger at apoapsis. 

Tab. [3] summarises the results of the complete study. 
The mean values of the eccentricity e, the error of the single 
value s e and the error of the mean value Se have been calcu- 
lated over time intervals in which the eccentricity remains 
stable, and are based on output data of a step- width of 0.002 
N-body units. The quantity etrans represents a mean value of 
the eccentricity in a period the total energy snaps of into the 
stage of constant energy loss. Although the the eccentricity 
still undergoes vigorous fluctuations at that time, this tran- 
sition mean value etrans does not evidently differ from the 
e. It is remarkable that eccentricity develops in a dynamical 
friction dominated period. The strength of the dynamical 
friction force along the trajectory of the black holes is cru- 
cial whether the individual motion determined by the initial 
velocities can be kept until the binary forms. 

The existence of drifts in eccentricity was reported in 
previous papers (see introduction). This behaviour can be 
connected to the effect of super-elastic scattering events. 
Although the effect of super-elastic scatterings on the ec- 
centricity during the formation of the binary may be only 
weak, the long term evolution might increase the eccentricity 
significantly. However, the question, whether the eccentric- 
ity is ultimately increased or decreased, cannot be solved 
here since the simulations includes both rising and declin- 
ing drifts, as seen for Wo — 6, uio = 0.0 with vo = v c and 
wo = \plv c (Fig|3j for example. If superelastic scattering 
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Figure 3. Evolution of the eccentricity of the IMBH binary for the complete survey. Starting from the top and from the left to the 
right, wc display in each panel only three models for readability. In each individual panel, and from the top to the bottom at t = 0, we 
show (in red, green and blue in the on-line version): A, B, C in the first panel, J', K, C in the second one, T>, £, T in the third one, M, 
Af, O in the fourth one, Q, H, I in the fifth one and V, Q in the last one 



should actually increase the eccentricity in long-term evolu- 
tion, another effect overlays the simulations in that period. 

We depict in Fig. [6] the cases T, C, 1, O, C and Q of 
Fig. [3] together in order to illustrate clearly the influence of 
rotation on the development of eccentricity. With a rotation 
parameter ojq — 0.6, the end value of eccentricity drops sig- 
nificantly for Wo — 3 as well as for Wo = 6. In the case of 
Wo = 3, u>o = 0.6, no effective transfer of angular momen- 
tum to the field stars takes place until t ~ 5, with Wo = 6, 
loo ~ 0.6 till t ~ 2, thus the trajectories of the IMBHs ex- 
perience a noticeable truncation in the co-rotating stellar 
system at the beginning, which finally results in low eccen- 
tricities. 



4.3 Inclination and angular momentum 
orientation 

The direction of the orbital angular momentum vector of 
the binary is specified by the quantities 



cosi = l-e z /l < i < ty (9) 

and 

{ fi = K = 0^fi<2^ (10) 

where e^ and e z are the unit vectors of the reference co- 
ordinate system and K the vector in the direction of the 
ascending node, K = e z x 1, which represents the definition 
of the inclination i and the longitude of the ascending node 
Q. following classical celestial mechanics. 

The left side of Fig. [7] shows the evolution of the incli- 
nation for Wo = 6 and u>o = 0.6 for different initial veloci- 
ties. The inclination, typically in all simulations undergoes 
comparatively strong changes during the dynamical friction 
dominated regime, and remains passably stable or slightly 
drifting during the hardening stage. Nevertheless, the total 
changes of the inclination angle considering the whole simu- 
lation are rather small. Sometimes a monotonic increase of 
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Figure 4. Projection in the X-Y plane of the trajectories of the two IMBHs in the model £ (different colours depict different IMBHs). 
Solid lines indicate the trajectories passed through in the time interval mentioned above each figure; the dotted lines show the previous 
orbit. Note the scaling of the axes in different figures; the two lower panels are a zoom in 



the inclination angle occurs during the dynamical friction 
dominated stage, and there were also simulations in which 
the inclination dropped before reaching stable behaviour. 
This can be seen in Tab. [4] where the maximum inclination 
imax and the mean value seen in the stable phase i of each 
simulation are listed: i max and i can differ considerably. 

In the right side of Fig. [7] the direction of the orbital 
angular momentum vector is illustrated using polar repre- 
sentation (i cos Q, i sin f2). The (0,0) coordinate corresponds 



to a rotation of the binary in the zy-reference plane. The ex- 
tensive ripples are the result of periodic motion before the 
binary is bound or when the binding is weak. With progress- 
ing evolution, the system concentrates in confined cloud-like 
areas. The measured changes of direction of the orbital an- 
gular momentum vecto r are consistent in order of mag nitude 
with the simulations of lMilosavlievic fc Merrittl l|2001l ), with 
the exception of the inclination maverick Wo — 6, coo — 0.0, 
Wo = 0.136u c . 



>, 
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Figure 5. the same as in Fig|4] but for model T 



5 COUNTER ROTATION 

Applying exactly the same initial model for the field stars of 
the corresponding model, the initial velocities of the IMBHs 
were set contrary the rotation of the stellar system for the 
models O and M ■ 

The distinct influence of dynamical friction during the 
the first time intervals is responsible for different results 
in eccentricity evolution compared to the co-rotating sim- 
ulation (Fig. [8}. In the case of vo = 0.136« c , for counter- 
rotation an extremely high eccentricity e = 0, 997 is reached 



(e = 0, 728 for co-rotating). In this scenario, the relative ve- 
locity between a black hole and the field stars is increased at 
the apoapsis as well as at the periapsis. Dynamical friction 
is very efficient at the apoapsis of the individual black hole 
motion in the unbound regime. Thus the black holes suffer a 
strong energy loss and fall steeper to the centre than in the 
co-rotating simulation. The resulting high eccentric motion 
is kept into the bound state. Applying an initial velocity 
Wo = Vc, a higher eccentricity e = 0.160 occurs compared to 
the co-rotating e = 0.039, but remains on a low level. 
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Figure 6. Influence of rotation on final eccentricities. Left panel: King Wo = 3, models T, C and T in green, red and blue, respectively. 
Right panel: King Wo = 6, models O, C and Q in green, red and blue, respectively. The colours indicate the different rotation parameters 
used in the simulations: u>o = 0.0 red, uio = 0.3 green, u>o = 0.6 blue 
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Table 3. Compiled data of eccentricity evolution for the the complete set of simulations. t sta b represents the time when the eccentricity 
remains passably constant, e is the average eccentricity measured over the time interval mentioned below, e en( j the value at the end of 
the simulation. Disregarding appearing drifts, errors were calculated based on an 0.002 step-width output considering the same time 
intervals; s e is the error of the single value and etrans the error of the mean value eccentricity, etrans is an average eccentricity over the 
epoch when the system turns to a constant energy loss rate 



6 BROWNIAN MOTION 

The centre of mass (CM) of a hardened binary is expected 
to to perform an irregular motion in the central region of the 
stellar system. This motion is often described by the concept 
of Brownian motion, as it is characterised by a friction force 
(dynamical friction) and a fluctuating force (as the result of 
scattering events and encounters of field stars). Applying en- 
ergy equipartition in thermodynamic equilibrium, the mean 
square velocity of the CM of the binary (i>cm) ls connected 
to the mean square velocity of the central field stars by 



(«cm) 



M 2 BH 

where M2BH is the sum of the black hole masses. 



(11) 



However, the irregular motion of the CM is to be dis- 
tinguished from the movement of a single massive particle 
since the binding energy of the binary changes due to (super- 
elastic) scattering events. The characteristics of the Brow- 
nian motion of a massive black h ole binary have been dis- 
cussed in detail bv lMerrittl (|200ll ). where (vq M ) is expected 
to be increased by a factor < 2, allowing for the higher re- 
coil velocities of a binary after super-elastic scattering of 
field stars and for the decreased dynamical friction force on 
the CM, since the trajectories of the field stars are randomly 
orientated in direction after such a process. 

The CM motion was investigated for a series of King pa- 
rameters. Fig(9]displays the CM movement during the whole 
simulation time for King potential Wo ~ 6, applying rota- 
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Figure 8. Left panel: Influence of rotation and counter-rotation on the evolution of the eccentricity of the IMBHs. The upper curve (at 
later times, green in the on-line version of the article) corresponds to model 0, for IMBHs initially set up to be counter-rotating; the 
lower curve (red colour) is the same model with co-rotation of the IMBHs. Right panel: Same-same for model Af 



tion parameters uo = 0.0 and uo = 0.6, both with an initial 
IMBH velocity vq = 0.136i> c . A color gradient plot of the 
CM trajectory has been used to give an impression of the 
temporal evolution. In the absence of rotation luq = 0.0, the 
characteristics of an rather irregular motion appears, while 
with rotation u>o = 0.6, a turning motion of the CM oc- 
curs. This motion seems to follow the rotation of the stellar 
system. 

Elevated values of the CM mean velocity were found 
in all simulations, quantified by the ratio (i>cM)/( w equ) with 
(^equ) = (h-*/A^2bh)o"o- The results of three simulations are 
shown in Tab(5] The postulated factor < 2 is manifestly 
exceeded even without rotation, but may be smaller if a a 
change of the central stellar velocity dispersion a is taken 
into account, while calculations here are based on the ini- 
tial value (To within the 1% Lagrangian radius of the stellar 
model. Nevertheless, a time evolution of the mean square 
velocity indicated that a constant value was not yet reached 



at the end of the simulations but may have grown larger if 
runs had been continued. 



7 GRAVITATIONAL WAVES: 

DETECTABILITY OF THE SYSTEMS WITH 
LISA 

7.1 Post evolution of the binary of IMBHs 

The advantage of direct-summation codes, accuracy, is at 
the price of performance. We have chosen N— body in or- 
der to investigate this problem but in order to analyse 
the ulterior evolution of the binary down to a GW fre- 
quency observable by LISA, we have to resort to alterna- 
tive schemes. If we were to integrate the binary system un- 
til the orbital period of the binary is within the range of 
observations of LISA, we would have to leave the simula- 
tions running for months. This is not desirable for obvious 
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Figure 9. Upper panel: Centre-of-mass trajectory of model C (u>o = 0.0). Lower panel: Same for model Q (uio = 0.6) The cross indicates 
the moment in which there is a constant energy loss rate (t ~ 7.6 for luq = 0.0 and t ~ 8.4 for ujq = 0.6) 
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Figure 10. Inspiral of the IMBH binary of Models C,F,C and O from the top to the bottom and from the left to the right. We show 
in this figure the evolution of the binary in the eccentricity— semi-major axis plane. The zigzag curves (red in the on-line version) depict 
the results of the N— body simulations. Initially the IMBHs do not constitute a bound system and therefore the eccentricity of the 
system is ill-defined, as explained at the beginning of Section l|4,ljl . This is causing the initial values of e become > 1. We then make a 
semi-analytical expansion of the evolution starting from the last point in the numerical simulations (dot-dashed curve, magenta in the 
on-line version of the paper) by taking into accound the Physical propierties of the system at that moment: semi-major axis, eccentricity, 
velocity dispersion of the stellar system and stellar density, as shown in Table [6] Thus, we evolve the binaries until they enter the LISA 
bandwidth, which we depict as a li ghtly shaded area (orange in the on-line version). The black solid curves correspond to the trajectories 
due only to the emission of GWs |Peterslll964l) and we additionally show the corresponding inspiral *gw f° r 10 10 yrs, 10 9 yrs etc. As 
shown in ( Amar o-Seoane & Freitag 2006), one recovers partially the N— body results with the semi-analytical approach if one starts at 
a previous point in the curve corresponding to the numerical simulation. The black-shaded region on the right corresponds to the last 
stable circular orbit. 



reasons. Instead, we recur to a semi-analytical method to 
evolve the orbital elements of the binary taking into ac- 
count the dyn amics and the GW emission of t he system, as 
introduced in lAmaro-Seoane fc Freitad l)2006l ): we stop the 
direct-summation calculation after the initial strong fluctu- 



ating phase; when the eccentricity achieves a steady value. 
In order to locate in the LISA bandwidth the position of the 
binary, we employ the results of the direct-summation simu- 
lation and extend them with a semi-analytical method. The 
dynamics will, in general, tend to increase the eccentricity, 
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Table 4. Inclination characteristics regarding all simulations, 
imax denotes the maximum value at the mentioned time unit, i 
is the mean value over the mentioned time interval, which ranges 
from the time unit when a constant energy loss rate occurs until 
the end of the simulation. The mean values are calculated despite 
some runs showing drift, in which case they are indicated with a 
"d". 
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Table 5. Mean square velocities for the CM of the binary (»g M ) 
compared to the equilibrium mean square velocity (v^ qu ) of a 
single particle, for some simulations. The (vq m ) were calculated 
over the blue epochs of section $9$ and over [11.60;30.58] in the 
case Wo = 6, u>o = 0.3, vq = v c . 



whilst the emission of GWs circularies the orbit. These two 
processes are competitive. The basic idea to further evolve 
the binary is to split the evolution of both the semi-major 
axis and the eccentricity in two contributions, one driven by 
the dynamical interactions with stars (subscript dyn) and 
another by the emission of GWs (subscript GW), 



da 
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/ dyn 



de 
dt 



de\ 
It) 



+ 



dyn 
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GW 

(12) 



The G W terms are as given by the approach used in lPetera 
(1 19641 ). As fo r the dynamica l part, we resort to the scheme 
described in iQuinlanl l|l996j). For more details about th is 
approach, see Section 3 of lAmaro-Seoane &: Freitad (| 20061 ). 

In Fig. (|10p we depict the results of this approach for 
the models C,J-,C and O. We have chosen them for being 
the most interesting ones from a dynamical point of view, 
since they have the largest eccentricities by the moment of 
entrance in the observatory's sensitivity window. Also, from 
the standpoint of detectability, these are the most challeng- 



ing ones due to the same reason. We can see that in all four 
models the binary of IMBHs enters the bandwith with a 
residual eccentricity which, even if rather mild, is not negli- 
gible, from the point of view of detection. 

7.2 Event rates 



iFregeau et alj l|2006j ) estimated the detection of binary 
IMBHs by LISA for the single cluster channel. They assumed 
that any cluster undergoing a collisional runaway, such as 
those found in the Monte Carlo numerical simulations of 
iGiirkan et al.1 |2006), form a two very massive stars. These 
evolve separately and eventually may collapse and build two 
IMBHs separately, so that the IMBHs do not coalesce in the 
pro cess of their formation, bu t are b orn independently. 

lAmaro-Seoane k, Freitael (|2006l ) calculated the event 
ra te for formation of bi naries of IMBHs based on the results 
of lFregeau et all l|2006l ) in the scenario of two colliding clus- 
ters, the double cluster channel (see their Section 4, also for 
a more detailed explanation of the following events). Their 
work assumed that the IMBHs were already present at the 
centres of the two clusters undergoing the crash. When the 
two clusters merge, the IMBHs are drawn to the centre due 
to dynamical friction and constitute a binary which even- 
tually c oalesce. The difference in the c alcul ation of event 
rates of lAmaro-Seoane fc Freitad (|2006l ) and lFregeau et all 
(2006) is the number of IMBHs formed per cluster and the 
requirement for the host clusters to merge. 

In both estimations, the probability that a cluster 
evolves to the runaway phase was set to 0. 1 as an illus- 
trative case, though it can be as high as 0.5 (|Freitag et al.l 
2006b). Both works assumed that a runaway always leads to 
t he formation of an IMBH. We re fer t he reader to Se c tion 4 
of lAmaro-Seoane fc Freitael (|200rj ) and lFregeau et al.l (|2006l ) 
for a detailed explanation and expositi on of the uncertainties 
in the calculations. To summarise, the IFregeau et al.l (|2006l ) 
results for the LISA detector are 



rprcg|opt G [200, 250]yr _1 (13) 

r Fr0 g|pc S G [40, 50] yr" 1 . 

Where the subscript "Freg" stands for IFregeau et all (2006), 
the subscript "opt" for the optimistic estimation, assum- 
ing that the probability for a cluster evolving to the 
runaway stage is 0.5, and the subscript "pes" stands 
for pesimistic, which is the r esult of using 0.1 instead. 
lAmaro-Seoane fc Freitael l|2006j ) find the following results, 
where we use a nomenclature similar as above and set the 
probability for the host clusters to merge to 1 (these would 
decrease by a factor 10 if one was to use 0.1 instead, see 
discussion about the "UCDG channel" in their work) 



r AS F|o P t G [100, 125] yr" 1 
TasfIpos G [4, 5]yr _1 . 



(14) 



So that the contribution to the total number of binaries of 
IMBHs from both channels is 



Ttotlopt G [300, 375] yr" 
Ttot lpcs G [44, 55] yr _1 



(15) 
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Simulation 


Time 


Rc 


N 


<4 


2 




n 


m 




P 


p* 


atot* 


C 


36.0 


0.376 


7778 


0.255 


0.304 


0.929 


34878.9 


1.82* 10" 


5 


0.635 


0.544 


0.860 


T 


33.0 


0.376 


7627 


0.253 


0.261 


0.880 


34201.8 


1.82* 10" 


-5 


0.622 


0.534 


0.943 


£ 


25.0 


0.216 


4209 


0.476 


0.782 


1.428 


99739.3 


2.03* 10" 


-5 


2.025 


1.555 


1.575 


O 


12.0 


0.223 


4270 


0.383 


0.265 


0.956 


91828.0 


2.03* 10" 


-5 


1.864 


1.433 


0.967 



Table 6. All values are given in N— body units. In the table R c is the core radius, N the number of stars within the core radius, including 
the two IMBHs, a^ R the radial velocity dispersion (squared), cr|, the tangetial velocity dispersion (squared, ID), atot = y &\ + 2<r|, the 
velocity dispersion (3D, including IMBHs), n the particle density within the core radius (including IMBHs), rh the average particle mass 
within the core radius (including IMBHs), p the average mass density within the core radius (including IMBHs), p* the average mass 
density within the core radius (excluding IMBHs), atot* the velocity dispersion within the core radius (3D, excluding 2BH) 



Or, in the unlikely very pesimistic situation of having the 
host clusters to merge with a 0.1 probability, we have the fol- 
lowing "optimistic-pesimistic" and "pesimistic-pesimistic" 
results 



Ttctl^t G [210, 262.5] yr" 1 (16) 

rtatllS 6 [40.4, 50.5] yr" 1 

These results are encouraging enough that we address 
the parameter estimation of the sources. We describe the 
methods and results in the next sections. 

7.3 The gravitational waveform 

To investigate the detectability of these sources for LISA, 
we use the restricted post-Newtonian (PN) approximation 
for the GW, where we assume 2-PN corrections to both the 
conservative and adiabatic dynamics of the system, but con- 
serve the amplitude at the dominant Newtonian order. With 
this in mind, the waveform polarisations for non-spinning 
eccentric bi naries at the detecto r are given by (in units of 
G = c = 1) jDamour et all j2004 )) 



h+(t) 



h x (t) 



Mr\ 



1 + COS I 

+ 



D L 

+ 2rr if sin 2ip 
2Mri 

= — COS L 

D L 

— 2rr if cos 2<y3 



cos 2ip I — f 2 + r 2 ip 2 + 



M 



■ 2 2 • 2 , M . j , _ 

-r — r if H J sm t f(17) 



.2 . 2 -2 , M 

-r + r if H sm lip 

r 



(18) 



where cos i — L n. Here, L is the direction of the binary's or- 
bital angular momentum and n is the direction from the ob- 
server to the source (such that the GWs propagate in the — n 
direction) . M is the total mass of the system, r\ — Mi M2 /M 2 
is the symmetric mass ratio and Dl is the luminosity dis- 
tance to the source . The components (r, <f) denote the or- 
bital radius and phase of the system (also referred to as 
the true a nomaly), and are sc hematically described by the 
equations l|Hinder et all (120081 )1 



r/M 
M<j> 
I 

Mi 



ronsix 1 + rip N + r 2P NX + 0(x 2 ) (19) 

+ 0(x 9/2 ) (20) 
(21) 

Mn = x 3/2 + /i PN x 5/2 + ! 2 pn/ /2 + 0(x 9/2 )(,22) 



1 3/2 . : 5/2,1 7/2 

<P0PNX + IplPNX + <p2PNX 

hpN + hpNX 2 + 0(x 3 ) 



where / is the mean anomaly and n — 2n/P is the mean mo- 
tion, where P is the orbital period defined as the time to go 
from pericenter to pericenter. Due to precession effects, this 
is different from the time taken to go from if to cp + 2tt. The 
dot represents the time derivative, e.g. if = dtp/dt. We also 
define the invariant PN coefficient x = (Ma)) 2 ' 3 , where the 
angular frequency is defined as u = (2-7T + Aip) /P and Aif is 
the precession angle of the pericenter per period. We should 
note that the coefficients in the equations presented above 
and below are in general functions of the instantaneous ec- 
centricity e and the eccentric anomaly u. In the limit of the 
eccentricity e 0, we reclaim the circular orbit case where 
u = (p. The adiabatic evolution of x and e are given by 



Mx = iopN^ 5 + iipN^ 6 + K2PN2: 7 + 0(x 
Me = eopr-jz 4 + eipNa; 5 + e2PNa; 6 + 0(x 



15/2n 



13/2s 



(23) 
(24) 



where we point the reader to the appendix of lllinder et all 
(2008) for the complete description of the above equations. 
The waveforms are generated as follows : after we have 
evolved the above equations for x and e, we numerically 
integrate Equation [221 to obtain l(t). We then use this value 
to solve the post-Newtonian Kepler's Equation 1211 which is 
a transcendental equation in u(t). Once we have u(t),e(t) 
and x(t), we can then calculate r(t),r(t),<p(t),ip(t), where 
we calculate the integral of <fi(t) and the derivative r(t) nu- 
merically. 

In this work, to fully describe the GW po- 
larisations we use the following parameter set : 
A = {ln(M c ),ln(^),ln(Dz,),ln(ao),cos6', 0,eo,cost, V 1 }, 
where M c = Mrj 3/l5 is the chirp-mass, ji = Mr\ is the 
reduced mass, Dl is the luminosity distance, ao is the 
initial semi-major axis of the orbit, (9, <j)) are the co- latitude 
and longitudinal position of the source in the sky, eo is the 
initial eccentricity, t is the inclination of the orbital plane 
and ip is the GW polarisation angle. 



7.4 Detector response and parameter error 
estimation 

LISA can be thought of as pair of co-located detectors ro- 
tated with respect to each other by an angle of 45 degree s. 
In the Low Frequency Approximation (LFA) (|Cutler| [l998). 
we can write the individual detector responses as 



hi(t) = h+ (mwHt) + h x mm* (t), 



(25) 
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where £(t) is the phase shifted time parameter defined by 

£(t)=t-R (B sin 6 cos (a(t)-c/>). (26) 

Here 7?® corresponds to one AU and a(t) = 2nf m t + k, 
where the LISA modulation frequency is f m = 1 / 'year and 
k is the initial ecliptic longitude of the guiding center of the 
LISA constellation. The functions F + ' x (t; 9, <j>, tjj) are the 
beam pattern functions of the detector (|Cornish fc Rubbol 

<Eqo3)). 

Given sources h(t) and g(t) we can define a noise 
weighted inner product 

df 



(A Iff) 



Sn(f) 



Hf)g*(f) + h*(f)~g(f). (27) 



where h(f) is the Fourier transform of the signal, an asterisk 
denotes a complex conjugate term and S n (f) is the one-sided 
noise spectral density of the detector. For this study, we use 
a noise curve given by 

s„(/) = sr tr (/) + sr f (/) (28) 

where the instrumental noise S n (f) is given by 



qui 



(29) 



I -f COS" ( j- 



(2tt10~ 



(2^/) 4 



(2tt/)« 



jCornish fc Rubbdl2003f ) . Here L 
length for LISA, ST S (/) = 4 x 10" 



5 x 10 km is the arm- 
2 m 2 /Hz and S„ c (f) = 



9x10 m /s /Hz are the position and acceleration noises 
respectively. The quantity /» = 1/ (2irL) is the mean trans- 
fer frequency for the LISA arm. Note that we have also in- 
cluded a reddened noise term which steepens the noise curve 
between 10 -4 and 10 -5 Hz. To model the galactic or con- 
fusion noise we use the following confusion noise estimate 
derived from a Nelem ans, Yungelson, Zwar t (NYZ) galac- 
tic fo reground model dNelemans et alll2004l ; iTimpano et all 
l200d ) 

-44.62 r-2.3 



10" 
10" 
10" 
10" 



7" 



62.8 



10~ 4 < / < 10" J 
10~ 3 < / < 10" 2 ' 7 

10" 2 ' 7 < / < io~ 2 - 4 

10" ~ 



(30) 



-89.68 f -20 10- 2 - 4 < / ^ 10- 2 

where the confusion noise has units of Hz" 1 . 

Using the noise weighted inner product, we can define 
the optimal signal to noise ratio (SNR) by 

df 



{h\h) c 



\h\})\ 



(31) 



/opt 'Jo S n (f) 
We can also define the Fisher information matrix (FIM) by 

T aS = {d a h\dph) . (32) 

where the theoretical standard deviation error estimate in 
parameter recovery is given as 



- \/ (r m ) (33) 

The derivatives of the waveforms appearing in the 
FIM are generated num erically. We refer the reader to 
l|Porter fc Cornish! [2008) for the intricacies in the numeri- 
cally calculation of the FIM. 



7.5 Sampling the parameter space 

Figure 10 provides the semi-major axes, eccentricities and 
orbital periods for the models C,J-,£ and O, assuming an 
equal mass IMBH binary with individual rest-frame masses 
of Mi = 440 Mq . To sample the parameter space we ran 
a 1000 iteration Monte Carlo simulation over the param- 
eters {ao, eo, 0, <t>, l, while keeping the luminosity dis- 
tance and redshifted mass parameters constant. For the 
angular parameters we assume (cos 8, cost) G [—1,1] and 
{(j>, ip) € [0, 2ir] and we then draw uniformly from these 
ranges. 

As the IMBHs in the study are quite low mass as com- 
pared to SMBHs, we need to ensure that the sources are 
detectable. Given the masses of the systems in question, the 
GW frequency at the last stable circular orbit is between 
4-6 Hz, which is well outside the frequency range of LISA. 
Therefore, we placed all sources at a common distance of 
Dl=100 Mpc and required a SNR greater than 5. At this 
distance, the parameter values observed by LISA are the red- 
shifted rather than rest frame values. To account for this, the 
measured total mass is M(z) — (1 + z)M and the measured 
GW frequency of the waveform is / gw (z) = /gw/(l + z). In 
this study we use the following relation between redshift, z, 
and luminosity distance, Dl ■ 



Dr 



c(l + z) 
Ho 



dz 



■si* 



+ n A \ 



where we assume WM AP values of values of (Qr , Om , ) = 
(4.9 x 10" 5 , 0.27, 0.73) and a Hubble's constant of #0=71 
km/s/Mpc. 

We also decided to enforce a maximum possible GW 
frequency of 3 mHz to e nsure the fidelity of the LFA 
l|Cornish fc Rubbol 120031 ; IShapiro Key fc Cornish! 120091 ). 
Thus, using the information in Figure 11, we evolved the 
sources from first entering the LISA bandwidth to the point 
where they achieved the required SNR threshold. For this 
we assumed a 3 year mission lifetime for LISA. 

Using the above constraints it was possible to find min- 
imum and maximum values of (oo,eo) which satisfied both 
the SNR and maximum frequency constraints. Furthermore 
to sample this parameter sub-space we found it was possible 
to relate the eccentricity and semi-major axis for the four 
models according to a quadratic law, i.e. 



-co + Ci ao + C2 a , 



(35) 



where we provide the coefficients d , the maximum and mini- 
mum values of both ao and eo in Table[7] Thus, by uniformly 
sampling ao, we also have a corresponding sample in eo. We 
can see from Table [7]that although the binaries have appre- 
ciable eccentricities when they first enter the LISA band, i.e 
e ~ 0.1 — 0.15, by time the systems become observable the 
eccentricity has dropped to e ~ 0.012 — 0.019. 

While the Monte Carlo is carried out using the sky co- 
ordinates of the source, a more interesting quantity to quote 
is LISA's angular resolution for each source. We define the 
angular resolution as 



(36) 



where 
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Model 


CO 




CI 


C2 


ag 11 " x 10 ~ 8 pc 


ag™ x 10~ s pc 


min 

e o 




e 




C 


1.1366 x 10" 


-3 


9.6414 x 10 4 


9.4069 x 10 11 


3.6123 


10.669 


3.5924 x 10" 


-3 


1.9942 x 10" 


-2 


T 


6.9767 x 10" 


-4 


6.3737 x 10 4 


6.9846 x 10 11 


3.6107 


9.8673 


2.5227 x 10" 


-3 


1.2389 x 10" 


-2 


£ 


1.1591 x 10" 


-3 


1.0606 x 10 5 


1.1698 x 10 12 


3.7049 


9.3739 


4.3852 x 10" 


-3 


1.9056 x 10- 


-2 


O 


8.6073 x 10" 


-4 


8.4893 x 10 4 


1.0181 x 10 12 


3.8524 


9.0204 


3.9227 x 10" 


-3 


1.5081 x 10" 


-2 



Table 7. Quadratic power law coefficients for the four models C ,T ', C,0 , as well as minimum and maximum ranges for the initial 
semi-major axes ag and eccentricity eg. 



E a/3 = ^AA a AA^ = (ivr 1 . (37) 

We can now define the quantities appearing in the angu- 
lar resolution as E 9e = (A cos A cos 6), E^ = (A(j>A(j)} 
and E e0 = (A cos 6 Aai) . The angular resolution has units of 
steradians. 

7.6 Results of the Monte Carlo 

In Figure [11] we plot the recovered SNRs for the four models 
of IMBH inspiral. The models C, T , C and O are represented 
top to bottom, and left to right. We can see that while it is 
possible to have strong sources, with SNRs 300, the vast 
majority of samples returned more modest, but detectable 
SNRs in the range of 5 to 50. This means that IMBH inspi- 
rals should be observable with the LISA detector. 

Due to the lower mass ranges of these binaries, for sys- 
tems placed at 100 Mpcs, the sources start to become visible 
in the detector at GW frequencies of 5 x 10 -4 Hz and higher. 
However, as they are still very widely separated at this point, 
there is very little evolution in frequency or eccentricity over 
a three year period. In Figure Q2] we plot the initial and fi- 
nal eccentricity distributions for the four models. We can see 
that there is very little circularisation of the binaries during 
the observation period, which means that systems entering 
the observable LISA band with eccentricities of ~ 0.02 will 
reach the end of the three years with almost the same ec- 
centricity. As a consequence, these sources should retain a 
measurable eccentricity throughout the LISA mission life- 
time. We discuss about the consequences of these results on 
lower-frequency Astrophysics and Data Analysis in the next 
section. 



8 CONCLUSIONS 

In this study we have carried out a dynamical study and a 
first step analysis of the detection of IMBH binary systems 
in rotating clusters. For the case of a rotating King model 
without rotation, th e results of the presented survey verify 
previo us outcomes bv lMakino et all l|l993l ). lHemsendorf et all 
( 2002]) for massive black hole binary evolution in Plummer 
models facing the development of the binding energy, the 
eccentricity and determined hardening constants. Analysing 
an extensive number of simulations, the main results from 
our study of the Dynamics of these systems can be described: 
(1) The final eccentricity is strongly dependent on the initial 
black hole velocities. (2) The eccentricity is dependent on 
the rotation parameter of the model. (3) Determined hard- 
ening rates in the same range of previous direct iV— body 



simulations of comparable particle numbers. (4) Only weak 
changes in the inclination and in the orientation of the an- 
gular momentum vector dir ection have been observed, con - 
sistent with simulations by iMilosavlievic fc Merrittl (|200lf ) . 
(5) Counter rotation simulations yield noticeable different 
results in eccentricity, in one case actually an extreme large 
value e = 0.997. (6) Brownian motion of the centre of mass of 
the binary is influenced by the rotation of the stellar system. 
All simulations indicate that the orbital parameters eccen- 
tricity and inclination develop to passably constant values 
in the non- or only weak bound state, determined by initial 
conditions and the influence of dynamical friction. 

In order to understand the impact of these sources in 
lower-frequency GW Astrophysics, we have extended the di- 
rect N— body simulations with a simplified semi-analytical 
model. Whilst this approach is a "kludge" and can only be 
envisaged as an approximation, the integration of the sys- 
tem down to LISA's window is out of question because this 
would require months of CPU calculation and, on the top 
of that, the numerical error would accumulate, so that the 
results would not be as robust as what one can expect from 
direct-summation schemes. 

We choose the systems yielding a larger eccentricity in 
the dynamical simulations because these are the most ap- 
pealing cases in the sense that their detection will be very 
challenging. Also, information about the previous dynamical 
story of the system is encoded in the radiation in the form 
of a nonnegiigible eccentricity. 

The results presented above show that LISA should 
have no problems in identifying the existence of IMBH bi- 
naries. Such events are important for LISA data analysis 
as they are a previously unconsidered source in terms of 
dectability and parameter extraction. Our simulations also 
suggest that they will spend their lifetime in the detector 
with a measurable eccentricity. In this work, we have looked 
at a particular case study where the masses and the luminos- 
ity distance of the sources were fixed, and a Monte Carlo ran- 
domisation carried out over the other response parameters. 
We demonstrated that we will be able to accurately measure 
the masses and sky resolution of such sources. While the ec- 
centricity is weak when the source becomes observable in 
the detector, it should still be possible to carry out a precise 
measurement of the initial eccentricity of the source. 

For this we used the LFA response for the LISA detec- 
tor. This limited the sources we investigated to a maximum 
GW frequency of 3 mHz to ensure that the LFA was still 
valid. As these are also quite high frequency sources, they 
have a long generation time, which puts a time constraint 
of the size of the Monte Carlo that can be carried out. Fi- 
nally, the waveforms used in this work represent eccentric 
non-spinning binaries. 
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Figure 11. The optimal SNRs for the inspiral of IMBH binary of Models C, T ', C and O from the top to the bottom and from the left to 
the right. While it is possible to have very strong sources with SNRs in the range of 300-400, the vast majority of samples return SNRs 
of between 5 and 50 for all four models 
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17.86 
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Table 8. Median values of the parameter estimation errors and SNRs for the four models C, T ', C,0 . Note that the units of AH is 
steradians 



As well as detectability, the extraction of the system 
parameters is also of great importance in GW Astrophysics. 
Using the FIM, we obtained the error predictions for the 
important system parameters. As the error predictions are 
a function of the position of the source in the sky, plus the 
orientation of the system with respect to the LISA constel- 
lation, the Monte Carlo simulations produced error distribu- 
tions with large tails. Because of this fact, we have decided 



to quote the median errors for the relevant parameters. In 
Table [8] we present the median errors for the parameters 
(M c , n, Dl, Af2, cos l, eo). We can see that for all models the 
fractional errors in the estimation of the chirp-mass and re- 
duced mass are of the order of 10 -4 and 10 -3 . While there 
is not much frequency evolution for these sources, the fact 
that they appear in the detector at frequencies on the order 
of mHz means that we can obtain errors in the luminosity 
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Figure 12. The initial and final eccentricities for the inspiral of IMBH binary of Models C,T,C and O from the top to the bottom and 
from the left to the right, assuming a 3 year observation period for LISA. As the binary components are still quite widely separated, 
there is little circularisation during the observation period, and the binaries thus retain a measurable eccentricity. Note the different 
eccentricity scales in the different cells 



distance of the order of 10 _1 . We see a similar order of er- 
ror in the estimation of cos t which is in general a difficult 
quantity to measure using electromagnetic information. 

It is also quite remarkable to see that angular resolu- 
tion of the IMBH inspirals is very good, with median errors 
on the order of 10 steradians, corresponding to an error 
box on the sky of about 3 square degrees. This level of accu- 
racy would place an inspiralling IMBH firmly in the field of 
view of a future detector such as the Large Synoptic Survey 
Telescope (LSST). Finally, we can also see, again from the 
fact that these sources are emitting GWs at frequencies on 
the order of mHz, the fractional errors in the estimation of 
initial eccentricity is on the order of 10 -7 . 

While we have shown that these IMBH binaries are de- 
tectable, there are a number of ways in which the analy- 
sis can be improved in the future. Firstly, a more repre- 
sentative study would also have randomised the individual 
masses of the binaries, as well as their luminosity distance. 
This would allow us to give a more concrete statement on 
detection and parameter estimation with the LISA detec- 



tor. Secondly, we restricted the maximum GW frequency 
of the binary to 3 mHz to ensure a valid approximation 
to the LISA response. In the future, we could investigate 
higher frequency binaries by either using a Rigid Adiabatic 
Approximation dRubbo et all 120041 ) or full response to the 
LISA detector. However, we should point out that for the 
higher frequency binaries, the initial eccentricity drops off 
rapidly, and these binaries may now be essentially circular. 
A more realistic study would also include the use of more 
realistic waveforms which include spins and higher harmon- 
ics. However, work on such templates has yet to fully start 
in earnest. Finally, it would also be interesting to carry out 
a longer Monte Carlo, and assume different mission lifetimes 
to see how detectability changes over observation time. 
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